<?php

/**
 * @package JAMA
 *
 *    For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n
 *    unit lower triangular matrix L, an n-by-n upper triangular matrix U,
 *    and a permutation vector piv of length m so that A(piv,:) = L*U.
 *    If m < n, then L is m-by-m and U is m-by-n.
 *
 *    The LU decompostion with pivoting always exists, even if the matrix is
 *    singular, so the constructor will never fail. The primary use of the
 *    LU decomposition is in the solution of square systems of simultaneous
 *    linear equations. This will fail if isNonsingular() returns false.
 *
 * @author  Paul Meagher
 * @author  Bartosz Matosiuk
 * @author  Michael Bommarito
 * @version 1.1
 * @license PHP v3.0
 */
class PHPExcel_Shared_JAMA_LUDecomposition {

	const MatrixSingularException = "Can only perform operation on singular matrix.";
	const MatrixSquareException = "Mismatched Row dimension";

	/**
	 *    Decomposition storage
	 *
	 * @var array
	 */
	private $LU = array();

	/**
	 *    Row dimension.
	 *
	 * @var int
	 */
	private $m;

	/**
	 *    Column dimension.
	 *
	 * @var int
	 */
	private $n;

	/**
	 *    Pivot sign.
	 *
	 * @var int
	 */
	private $pivsign;

	/**
	 *    Internal storage of pivot vector.
	 *
	 * @var array
	 */
	private $piv = array();


	/**
	 *    LU Decomposition constructor.
	 *
	 * @param $A Rectangular matrix
	 *
	 * @return Structure to access L, U and piv.
	 */
	public function __construct($A) {
		if ($A instanceof PHPExcel_Shared_JAMA_Matrix)
		{
			// Use a "left-looking", dot-product, Crout/Doolittle algorithm.
			$this->LU = $A->getArray();
			$this->m = $A->getRowDimension();
			$this->n = $A->getColumnDimension();
			for ($i = 0; $i < $this->m; ++$i)
			{
				$this->piv[$i] = $i;
			}
			$this->pivsign = 1;
			$LUrowi = $LUcolj = array();

			// Outer loop.
			for ($j = 0; $j < $this->n; ++$j)
			{
				// Make a copy of the j-th column to localize references.
				for ($i = 0; $i < $this->m; ++$i)
				{
					$LUcolj[$i] = & $this->LU[$i][$j];
				}
				// Apply previous transformations.
				for ($i = 0; $i < $this->m; ++$i)
				{
					$LUrowi = $this->LU[$i];
					// Most of the time is spent in the following dot product.
					$kmax = min($i, $j);
					$s = 0.0;
					for ($k = 0; $k < $kmax; ++$k)
					{
						$s += $LUrowi[$k] * $LUcolj[$k];
					}
					$LUrowi[$j] = $LUcolj[$i] -= $s;
				}
				// Find pivot and exchange if necessary.
				$p = $j;
				for ($i = $j + 1; $i < $this->m; ++$i)
				{
					if (abs($LUcolj[$i]) > abs($LUcolj[$p]))
					{
						$p = $i;
					}
				}
				if ($p != $j)
				{
					for ($k = 0; $k < $this->n; ++$k)
					{
						$t = $this->LU[$p][$k];
						$this->LU[$p][$k] = $this->LU[$j][$k];
						$this->LU[$j][$k] = $t;
					}
					$k = $this->piv[$p];
					$this->piv[$p] = $this->piv[$j];
					$this->piv[$j] = $k;
					$this->pivsign = $this->pivsign * -1;
				}
				// Compute multipliers.
				if (($j < $this->m) && ($this->LU[$j][$j] != 0.0))
				{
					for ($i = $j + 1; $i < $this->m; ++$i)
					{
						$this->LU[$i][$j] /= $this->LU[$j][$j];
					}
				}
			}
		}
		else
		{
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
		}
	} //	function __construct()


	/**
	 *    Get lower triangular factor.
	 *
	 * @return array Lower triangular factor
	 */
	public function getL() {
		for ($i = 0; $i < $this->m; ++$i)
		{
			for ($j = 0; $j < $this->n; ++$j)
			{
				if ($i > $j)
				{
					$L[$i][$j] = $this->LU[$i][$j];
				}
				elseif ($i == $j)
				{
					$L[$i][$j] = 1.0;
				}
				else
				{
					$L[$i][$j] = 0.0;
				}
			}
		}

		return new PHPExcel_Shared_JAMA_Matrix($L);
	} //	function getL()


	/**
	 *    Get upper triangular factor.
	 *
	 * @return array Upper triangular factor
	 */
	public function getU() {
		for ($i = 0; $i < $this->n; ++$i)
		{
			for ($j = 0; $j < $this->n; ++$j)
			{
				if ($i <= $j)
				{
					$U[$i][$j] = $this->LU[$i][$j];
				}
				else
				{
					$U[$i][$j] = 0.0;
				}
			}
		}

		return new PHPExcel_Shared_JAMA_Matrix($U);
	} //	function getU()


	/**
	 *    Return pivot permutation vector.
	 *
	 * @return array Pivot vector
	 */
	public function getPivot() {
		return $this->piv;
	} //	function getPivot()


	/**
	 *    Alias for getPivot
	 *
	 * @see getPivot
	 */
	public function getDoublePivot() {
		return $this->getPivot();
	} //	function getDoublePivot()


	/**
	 *    Is the matrix nonsingular?
	 *
	 * @return true if U, and hence A, is nonsingular.
	 */
	public function isNonsingular() {
		for ($j = 0; $j < $this->n; ++$j)
		{
			if ($this->LU[$j][$j] == 0)
			{
				return false;
			}
		}

		return true;
	} //	function isNonsingular()


	/**
	 *    Count determinants
	 *
	 * @return array d matrix deterninat
	 */
	public function det() {
		if ($this->m == $this->n)
		{
			$d = $this->pivsign;
			for ($j = 0; $j < $this->n; ++$j)
			{
				$d *= $this->LU[$j][$j];
			}

			return $d;
		}
		else
		{
			throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
		}
	} //	function det()


	/**
	 *    Solve A*X = B
	 *
	 * @param  $B  A Matrix with as many rows as A and any number of columns.
	 *
	 * @return  X so that L*U*X = B(piv,:)
	 * @PHPExcel_Calculation_Exception  IllegalArgumentException Matrix row dimensions must agree.
	 * @PHPExcel_Calculation_Exception  RuntimeException  Matrix is singular.
	 */
	public function solve($B) {
		if ($B->getRowDimension() == $this->m)
		{
			if ($this->isNonsingular())
			{
				// Copy right hand side with pivoting
				$nx = $B->getColumnDimension();
				$X = $B->getMatrix($this->piv, 0, $nx - 1);
				// Solve L*Y = B(piv,:)
				for ($k = 0; $k < $this->n; ++$k)
				{
					for ($i = $k + 1; $i < $this->n; ++$i)
					{
						for ($j = 0; $j < $nx; ++$j)
						{
							$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
						}
					}
				}
				// Solve U*X = Y;
				for ($k = $this->n - 1; $k >= 0; --$k)
				{
					for ($j = 0; $j < $nx; ++$j)
					{
						$X->A[$k][$j] /= $this->LU[$k][$k];
					}
					for ($i = 0; $i < $k; ++$i)
					{
						for ($j = 0; $j < $nx; ++$j)
						{
							$X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k];
						}
					}
				}

				return $X;
			}
			else
			{
				throw new PHPExcel_Calculation_Exception(self::MatrixSingularException);
			}
		}
		else
		{
			throw new PHPExcel_Calculation_Exception(self::MatrixSquareException);
		}
	} //	function solve()

} //	class PHPExcel_Shared_JAMA_LUDecomposition
